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Abstract 

Qs , Engineering optimization is typically multiobjective and multidisciplinary with complex con- 

^SJ ' straints, and the solution of such complex problems requires efhcient optimization algorithms. 

Recently, Xin-She Yang proposed a bat-inspired algorithm for solving nonlinear, global optimi- 
sation problems. In this paper, we extend this algorithm to solve multiobjective optimisation 
problems. The proposed multiobjective bat algorithm (MOBA) is first validated against a sub- 
^^ ' set of test functions, and then applied to solve multiobjective design problems such as welded 

-^ , beam design. Simulation results suggest that the proposed algorithm works efficiently. 

^ . 
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^ ■ 1 Introduction 

Design optimisation in engineering often concerns multiple design objectives under complex, highly 
nonlinear constraints. Different objectives often conflict each other, and sometimes, truly optimal 
solutions do not exist, and some tradeoff and approximations are often needed. Further to this 
^^ I complexity, a design problem is subjected to various design constraints, limited by design codes or 

?H ■ standards, material properties and choice of available resources and costs (Deb, 2001; Farina et al., 

2004). Even for global optimisation problems with a single objective, if the design functions are 
highly nonlinear, global optimality is not easy to reach. Metaheuristic algorithms are very powerful 
in dealing with this kind of optimization, and there are many review articles and excellent textbooks 
(Coello, 1999; Deb, 2001; Isasi and Hernandez, 2004; Yang, 2008; Talbi, 2009; Yang, 2010c). 

In contrast with single objective optimization, multiobjective problems are much difficult and 
complex (Coello, 1999; Floudas et al., 1999; Gong et al., 2009; Yang and Koziel, 2010). Firstly no 
single unique solution is the best; instead, a set of non-dominated solutions should be found in order 
to get a good approximation to the true Pareto front. Secondly, even if an algorithm can find solution 
points on the Pareto front, there is no guarantee that multiple Pareto points will distribute along the 
front uniformly, often they do not. Thirdly, algorithms work well for single objective optimization 
usually do not directly work for multiobjective problems, unless under special circumstances such 
as combining multiobjectives into a single objective using some weighted sum methods. Substantial 
modifications are often needed. In addition to these difficulties, a further challenge is how to generate 
solutions with enough diversity so that new solutions can sample the search space efficiently (Talbi, 
2009; Erfani and Utyuzhnikov, 2011; Yang and Koziel, 2011). 
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Furthermore, real-world optimization problems always involve certain degree of uncertainty or 
noise. For example, materials properties for a design product may vary significantly, an optimal 
design should be robust enough to allow such inhomogeneity and also provides good choice for 
decision-makers or designers. Despite these challenges, multiobjective optimization has many pow- 
erful algorithms with many successful applications (Abbass and Sarker, 2002; Banks et al., 2008; 
Deb, 2001, Farina et al., 2004; Konak et al., 2006; Rangaiah, 2008; Marler and Arora, 2004). 

In addition, metaheuristic algorithms start to emerge as a major player for multiobjective global 
optimization, they often mimic the successful characteristics in nature, especially biological systems 
(Kennedy and Eberhart, 1995; Yang, 2005; Yang, 2010a; Yang, 2010b). Many new algorithms 
are emerging with many important applications (Kennedy and Eberhart, 1995; Luna et al., 2007; 
Osyczka and Kundu, 1995; Reyes-Sierra and Coello, 2006; Tabli, 2009; Cui and Cai, 2009; Yang, 
2010c; Zhang and Li, 2007; Yang and Deb, 2010b, Yang et al., 2011). For example, a new cuckoo 
search algorithm was developed by Xin-She Yang and Suash Deb (2009) and more detailed studies 
by the same authors (Yang and Deb, 2010a) suggested that it is very efficient for solving nonlinear 
engineering design problems. For a recent review of popular metaheuristics, please refer to Yang 
(2011). 

Recently, a new metaheuristic search algorithm, called bat algorithm (BA), has been developed 
by Xin-She Yang (2010a). Preliminary studies show that it is very promising and could outperform 
existing algorithms. In this paper, we will extend BA to solve multiobjective problems and formulate 
a multiobjective bat algorithm (MOB A). We will first validate it against a subset of multiobjective 
test functions. Then, we will apply it to solve design optimization problems in engineering, such as 
bi-objective beam design. Finally, we will discuss the unique features of the proposed algorithm as 
well as topics for further studies. 

2 Bat Behaviour and Bat Algorithm 

In order to extend the bat-inspired algorithm for single optimization to solve multiobjective problems, 
let us briefly review the basics of the bat algorithm for single objective optimization. Then, we will 
outline the basic ideas and steps of the proposed algorithm. 

2.1 Echolocation of Microbats 

Bats are fascinating animals. They are the only mammals with wings and they also have advanced 
capability of echolocation. It is estimated that there are about 996 different species which account 
for up to 20% of all mammal species (Altringham, 1996; Colin, 2000). Their size ranges from the 
tiny bumblebee bat (of about 1.5 to 2g) to the giant bats with wingspan of about 2 m and weight 
up to about 1 kg. Microbats typically have forearm length of about 2.2 to 11cm. Most bats uses 
echolocation to a certain degree; among all the species, microbats are a famous example as microbats 
use echolocation extensively while megabats do not (Richardson, 2008). 

Microbats use a type of sonar, called, echolocation, to detect prey, avoid obstacles, and locate 
their roosting crevices in the dark. These bats emit a very loud sound pulse and listen for the 
echo that bounces back from the surrounding objects. Their pulses vary in properties and can be 
correlated with their hunting strategies, depending on the species. Most bats use short, frequency- 
modulated signals to sweep through about an octave, while others more often use constant-frequency 
signals for echolocation. Their signal bandwidth varies depends on the species, and often increased 
by using more harmonics. 

Though each pulse only lasts a few thousandths of a second (up to about 8 to 10 ms), however, 
it has a constant frequency which is usually in the region of 25kHz to 150 kHz. The typical range of 
frequencies for most bat species are in the region between 25kHz and lOOkHz, though some species 
can emit higher frequencies up to 150 kHz. Each ultrasonic burst may last typically 5 to 20 ms, and 
microbats emit about 10 to 20 such sound bursts every second. When hunting for prey, the rate of 
pulse emission can be sped up to about 200 pulses per second when they fly near their prey. Such 
short sound bursts imply the fantastic ability of the signal processing power of bats. In fact, studies 
shows the integration time of the bat ear is typically about 300 to 400 /is. As the speed of sound 



in air is typically v = 340 m/s, the wavelength A of the ultrasonic sound bursts with a constant 
frequency / is given by A = v/ f, which is in the range of 2nim to 14nim for the typical frequency 
range from 25kHz to 150 kHz. Such wavelengths are in the same order of their prey sizes. 

Studies show that microbats use the time delay from the emission and detection of the echo, 
the time difference between their two ears, and the loudness variations of the echoes to build up 
three dimensional scenario of the surrounding. They can detect the distance and orientation of the 
target, the type of prey, and even the moving speed of the prey such as small insects (Altringham, 
1996). Obviously, some bats have good eyesight, and most bats also have very sensitive smell sense. 
In reality, they will use all the senses as a combination to maximize the efficient detection of prey 
and smooth navigation. However, here we are only interested in the echolocation and the associated 
behaviour. Such echolocation behaviour of microbats can be formulated in such a way that it can 
be associated with the objective function to be optimized, and this makes it possible to formulate 
new optimization algorithms. 

2.2 Bat Algorithm 

If we idealize some of the echolocation characteristics of microbats, we can develop various bat- 
inspired algorithms or bat algorithms. In the basic bat algorithm developed by Xin-She Yang 
(2010a), the following approximate or idealized rules were used. 

1. All bats use echolocation to sense distance, and they also 'know' the difference between 
food/prey and background barriers in some magical way; 

2. Bats fly randomly with velocity Vi at position Xi with a frequency /min, varying wavelength 
A and loudness Aq to search for prey. They can automatically adjust the wavelength (or 
frequency) of their emitted pulses and adjust the rate of pulse emission r € [0, 1], depending 
on the proximity of their target; 

3. Although the loudness can vary in many ways, we assume that the loudness varies from a large 
(positive) Aq to a minimum constant value A^i^- 

Another obvious simplification is that no ray tracing is used in estimating the time delay and 
three dimensional topography. Though this might be a good feature for the application in compu- 
tational geometry, however, we will not use this feature, as it is more computationally extensive in 
multidimensional cases. 

In addition to these simplified assumptions, we also use the following approximations, for sim- 
plicity. In general the frequency / in a range [/min,/max] corresponds to a range of wavelengths 
[Amin, Amax]- For example a frequency range of [20kHz, 500kHz] corresponds to a range of wave- 
lengths from 0.7mm to 17mm in reality. Obviously, we can choose the ranges freely to suit different 
applications. 

2.3 Bat Motion 

For the bats in simulations, we have to define the rules how their positions Xi and velocities Vi in a 
d-dimensional search space are updated. The new solutions a;* and velocities v^ at time step t are 
given by 

Ji — ,/min i v,/max /minjPi \^) 

vl+'^vl + ixl-x^)f,, (2) 

xl+'^xl+vl (3) 

where /3 G [0, 1] is a random vector drawn from a uniform distribution. Here x^ is the current global 
best location (solution) which is located after comparing all the solutions among all the n bats at 
each iteration t. As the product Xifi is the velocity increment, we can use fi (or A^ ) to adjust the 
velocity change while fixing the other factor Xi (or fi), depending on the type of the problem of 
interest. In our implementation, we will use /min — and /max = 0(1), depending on the domain 



size of the problem of interest. Initially, each bat is randomly assigned a frequency which is drawn 
uniformly from [/min,/max]- 

For the local search part, once a solution is selected among the current best solutions, a new 
solution for each bat is generated locally using random walk 

a^now = a^oid + e A*, (4) 

where e is a random number vector drawn from [—1, 1], while A* =<A\ > is the average loudness of 
all the bats at this time step. 

The update of the velocities and positions of bats have some similarity to the procedure in the 
standard particle swarm optimization, as fi essentially controls the pace and range of the movement 
of the swarming particles. To a degree, BA can be considered as a balanced combination of the 
standard particle swarm optimization and the intensive local search controlled by the loudness and 
pulse rate. 

2.4 Loudness and Pulse Emission 

Furthermore, the loudness Ai and the rate r^ of pulse emission have to be updated accordingly as the 
iterations proceed. As the loudness usually decreases once a bat has found its prey, while the rate 
of pulse emission increases, the loudness can be chosen as any value of convenience. For example, 
we can use ^o = 100 and Amin — 1- For simplicity, we can also use Aq — \ and Amin — 0, assuming 
^min — means that a bat has just found the prey and temporarily stop emitting any sound. Now 
we have 

A'+^^aAl r*=rO[l-exp(-7t)], (5) 

where a and 7 are constants. In fact, a is similar to the cooling factor of a cooling schedule in the 
simulated annealing (Kirkpatrick et al., 1983). For any < a < 1 and 7 > 0, we have 

A\ ^ 0, r\ ^ rO, as i ^ c^. (6) 

In the simplest case, we can use a = 7, and we have used a = 7 = 0.9 in our simulations. 

Preliminary studies by Yang (2010a) suggested that bat algorithm is very promising for solving 
nonlinear global optimization problems. Now we extend it to solve multiobjective optimization 
problems. 

3 Multiobjective Bat Algorithm 

Multiobjective optimization problems are more complicated than single objective optimization, and 
we have to find and/or approximate the optimality fronts. In addition, algorithms have to be 
modified to accommodate multiobjectives properly. 

3.1 Pareto Optimality 

A solution vector u = (?ii, ..,m„) e T, is said to dominate another vector v — (wi, ...,w„) if and 
only if Ui < Vi for Vz € {1, ..., n} and 3i G {1, ..., n} : Ui < Vi. In other words, no component of u is 
larger than the corresponding component of v, and at least one component is smaller. Similarly, we 
can define another dominance relationship ^ by 

u ^ V '^=^ u -< vy u = V. (7) 

It is worth pointing out that for maximization problems, the dominance can be defined by replacing 
-< with )~. Therefore, a point a;* G J^ is called a non-dominated solution if no solution can be found 
that dominates it (Coello, 1999). 

The Pareto front PF of a multiobjective can be defined as the set of non-dominated solutions 
so that 

PF = {scS^s' C,S:s' <s}, (8) 



Objective functions fi{x), ..., fx{x), x = {xi,...,xj)^ 
Initialize the bat population Xi {i — 1, 2, ..., n) and Vi 
for j = 1 to A^ (points on Pareto fronts) 

Generate K weights Wk > so that J2k=i ^fe = 1 
Form a single objective / = J2k=i ^kfk 
while {t <Max number of iterations) 

Generate new solutions and update by ([T]) to ([3]) 
if (rand > r^) 

Random walk around a selected best solution 
end if 

Generate a new solution by flying randomly 
if (rand < A, & f{x^) < /(a;*)) 
Accept the new solutions, 
and increase r^ & reduce Ai 
end if 

Rank the bats and find the current best x^, 
end while 

Record a;* as a non-dominated solution 
end 
Postprocess results and visualization 



Figure I: Multiobjective bat algorithm (MOB A). 

or in term of the Pareto optimal set in the search space 

PF* = {xeT\^x'eT:f{x')^fix)}, (9) 

where f = {fi, ..., /k)^ ■ To obtain a good approximation to Pareto front, a diverse range of solutions 
should be generated using efficient techniques (Gujarathi and Babu, 2009; Konak et al., 2006). 

3.2 MOBA Algorithm 

Based on these approximations and idealization, the basic steps of the multiobjective bat algorithm 
(MOBA) can be summarized as the pseudo code shown in Fig. [TJ 

For simplicity here, we use a weighted sum to combine all objectives fk into a single objective 

K K 

f = Y,wkfk, j^wfc^i. (10) 

fc=i fc=i 

As the weights are generated randomly from a uniform distribution, it is possible to vary the weights 
with sufficient diversity so that the Pareto front can be approximated correctly. 

In our simulations, we have carried out parametric studies, and we have used a — "f = 0.9 for all 
simulations. The choice of parameters requires some experimenting. Initially, each bat should have 
different values of loudness and pulse emission rate, and this can be achieved by randomization. 
For example, the initial loudness A'^ can typically be [1,2], while the initial emission rate r^ can be 
around zero, or any value r^ S [0, 1] if using ([5]). Their loudness and emission rates will be updated 
only if the new solutions are improved, which means that these bats are moving towards the optimal 
solution. 



4 Numerical Results 

4.1 Parametric Studies 

The proposed multiobjective bat algorithm (MOBA) is implemented in Matlab, and computing 
time is within a few seconds to less than a minute, depending on the problem of interest. We have 
tested it using a different range of parameters such as population size (n), loudness reduction a, 
and pulse reduction rate 7. By varying n^5, 10, 15, 20, 30, 50 to 50, 100, 150, 200, 250, 300, 400, 500, 
a = 0.5,0.6,0.7, 0.8,0.9,0.95,1 and 7 = 0.5,0.6,0.7, 0.8,0.9,0.95,1,0, we found that the best 
parameters for most applications are: n = 25 to 50, a = 0.7 to 0.9 and 7 = 0.7 to 0.9. 

The stopping criterion can be defined in many ways. We can either use a given tolerance or a 
fixed number of iterations. From the implementation point of view, a fixed number of iterations is 
not only easy to implement, but also suitable to compare the closeness of Pareto front of different 
functions. So we have set the fixed number iterations as 5000, which is sufficient for most problems. 
If necessary, we can also increase it to a larger number. 

In order to generate more optimal points on the Pareto front, we can do it in two ways: increase 
the population size n or run the program a few more times. Through simulations, we found that to 
increase of n typically leads to a longer computing time than to re-run the program a few times. This 
may be due to the fact that manipulations of large matrices or longer vectors usually take longer. 
Another possibility is that simple restart can increase the diversity of solutions than more intensive 
search for longer iterations. So to generate 200 points using a population size n = 50 requires to run 
the program 4 times, which is easily done within a few minutes. Therefore, in all our simulations, 
we will use the fixed parameters: n ~ 50, a = 7 = 0.9. 

4.2 Multiobjective Test Functions 

There are many different test functions for multobjective optimization (Schaffer, 1985; Zhang et al., 
2003; Zhang et al, 2009; Zitzler and Thiele, 1999; Zitzler et al., 2000), but a subset of a few widely 
used functions provides a wide range of diverse properties in terms Pareto front and Pareto optimal 
set. To validate the proposed MOBA, we have selected a subset of these functions with convex, 
non-convex and discontinuous Pareto fronts. We also include functions with more complex Pareto 
sets. To be more specific in this paper, we have tested the following four functions: 
ZDTl function with a convex front (Zitzler and Thiele 1999; Zitzler et al. 2000) 



/i(x)=xi, J2{x) =.9(1 - \J !\lg), 

.9 = 1+ ^r^' ' 2;, G[0,1], z = l,...,30, (11) 

where d is the number of dimensions. The Pareto-optimality is reached when 5 = 1. 
ZDT2 function with a non-convex front 

/i(x)=xi, /2(x)-.g(l-:^)2, 

9 

ZDT3 function with a discontinuous front 

/i(x)=xi, f2{x) ^ g 1- J— - — sin(107r/i; 
L \j g 9 

where g in functions ZDT2 and ZDT3 is the same as in function ZDTl. In the ZDT3 function, /i 
varies from to 0.852 and /2 from —0.773 to 1. 

LZ4 function (Li and Zhang, 2009; Zhang and Li, 2007) 

2 ,-^ 



Table 1: Summary of results. 
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ZDTl 


3.7E-4 


4.5E-17 


ZDT2 


2.4E-4 


3.2E-19 


ZDT3 


5.2E-5 


1.7E-15 


LZ4 


2.9E-4 


1.2E-16 



/2 = 1 - X? 



Jt-'2 



/l(Uj), 



(12) 



where Ji = {j\j is odd and 2 < j < d} and J2 = {j\j is even and 2 < j < d}. 

Uj = Xj — sin(67rxi + — r): 



and 



XI € [0,1], x,& [-2,2], J =2,..., d, 



1 



»2|i;| 



This function has a Pareto front /2 = 1 — /f for < /i < 1 with a Pareto set 



?7r 
Cj = sin(67rxi + — ), j = 2, 3, ..., d, xi e [0, 1] 



(13) 



(14) 



After generating 200 Pareto points by MOBA, the Pareto front generated by MOBA is compared 
with the true front /2 = 1 — -y/Ti of ZDTl (sec Fig. [2]). In all the rest of the figures, the vertical 
axis is for /2 while the horizontal axis is for /i. 

Let us define the distance or error between the estimate Pareto front PF'^ to its correspond true 
front PF* as 



N 



E 



f 



(15) 



|p^e _ ppt||2 ^ J2{PFj - PFtf, 
J = l 

where N is the number of points. The convergence property can be viewed by following the iterations. 
Figs. [3] and [5] show the exponential-like decrease of Ef as the iterations proceed. The least-square 
distance from the estimated front to the true front of ZDTl for the first 1000 iterations (Fig. ^ and 
the logarithmic scale for 5000 iterations (Fig. S]). 

We can see clearly that our MOBA algorithm indeed converges almost exponentially. The results 
for all the functions are summarized in Table 1. We can see that exponential convergence can be 
achieved in all cases. 



5 Engineering Optimization 

Design optimization, especially design of structures, has many applications in engineering and in- 
dustry. As a result, there are many different benchmarks with detailed studies in the literature 
(Pham and Ghanbarzadeh, 2007; Ray and Liew, 2002; Rangaiah, 2008). Among the widely used 
benchmarks, the welded beam design is a well-known design problem. In the rest of this paper, we 
will solve this design benchmark using MOBA. 

Multiobjective design of a welded beam is a classical benchmark which has been solved by many 
researchers (Deb, 1999; Gong et al., 2009; Ray and Liew, 2002). The problem has four design 
variables: the width w and length L of the welded area, the depth d and thickness h of the main 
beam. The objectives are to minimize both the overall fabrication cost and the end deflection i5. 
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Figure 2: Estimated front and true front for ZDTl. 
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Figure 3: Convergence of the MOBA. 
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Figure 4: Exponential convergence of the MOBA. 



The problem can be written as 

minimise fi{x) = 1.10471w'^L + 0.04811rf/i(14.0 + L), 

minimize /2 — S, 
subject to 



where 



9i{x) 


^w- h<0, 




52(33) 


= S{x)-0.25<0, 




93(x) 


= t{x)- 13,600 <0, 




diix) 


= (7{x) -30, 000 < 0, 




55(33) 


:= 0.10471w2 + 0.04811M(14 + L) - 


-5.0< 0, 


96{x) 


= 0.125 -w < 0, 




57(33) 


= 6000 - Fix) < 0, 





a{x) 



504,000 



Q = 6000(14+1), 



X 65,856 a 

= T^TTT^fcn-^, P 



(w+df 



65,856 
30,000/id^ ' 



QD 
J ' 



6000 
\/2ji'L' 



r{x) = Ja^ + ^+l3' 



P = 0.61423 X 10" -^(1- 



dh^ dy/30/A8 



(16) 



(17) 



(18) 



(19) 



6 '" 28 ' 

The simple limits or bounds are 0.1 < L,d < 10 and 0.125 < w,h < 2.0. 

By using the MOBA, we have solved this design problem. The approximate Pareto front gener- 
ated by the 50 non-dominated solutions after 1000 iterations are shown in Fig. [SJ This is consistent 
with the results obtained by others (Ray and Liew, 2002; Pham and Ghanbarzadeh, 2007). In 
addition, the results are more smooth with fewer iterations. 

The simulations for these benchmarks and functions suggest that MOBA is a very efficient 
algorithm for multiobjective optimization. It can deal with highly nonlinear problems with complex 
constraints and diverse Pareto optimal sets. 
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Figure 5: Parcto front for the bi-objective beam design. 

6 Conclusions 

Multiobjective optimization problems are typically very difficult to solve. In this paper, we have 
successfully formulated a new algorithm for multiobjective optimization, namely, multiobjective bat 
algorithm, based on the recently developed bat algorithm. The proposed MOBA has been tested 
against a subset of well-chosen test functions, and then been applied to solve design optimization 
benchmarks in structural engineering. Results suggest that MOBA is an efficient multiojective 
optimizer. 

Additional tests and comparison of the proposed are highly needed. In the future work, we will 
focus on the parametric studies for a wider range of test problems, including discrete and mixed 
type of optimization problems. We will try to test the diversity of the Pareto front it can generate 
so as to identify the ways to improve this algorithm to suit a diverse range of problems. There are a 
few efficient techniques to generate diverse Pareto fronts (Erfani and Utyuzhnikov 2011), and some 
combination with these techniques may improve MOBA even further. 

Further research can also emphasize the performance comparison of this algorithm with other 
popular methods for multiobjective optimization. In addition, hybridization with other algorithms 
may also prove to be fruitful. 
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